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In order to determine the wind field in the lower atmosphere over 
the entire globe one common technique is to solve the "linear balance 
equation" , 

2 

fy -h . Vf = g($) (1) 

where ^ is the stress function whose gradient gives the wind field, f is 
the Coriolis function 

f = 2J2 sin <i> , (2) 

g($) is a known function of the (known) geopotential ^ ±s the earth’s 

angular velocity and (j) is the latitude (polar) angle measured from the 
equator. Equation Cl) is to be solved on the surface of a sphere of radius 
a and so is written 



2n . A 1 

^2 cos (}) 6(j) 



" (cos ♦ -JJ) + ^ 

cos (j) d A ) 



1 S 2^^ cos j) $4^ 



(5(}) 



g(^) . (3) 



In practical applications (3) is solved numerically, by finite difference 

methods, and the vanishing coefficient of the first term for small (}) causes 

real numerical difficulties, (The question of whether (1) is an appropriate 

equation to solve near the equator (for small <i>) is another argument to 

which we return later,) One conventional scheme to avoid these problems 

is to replace ^ constant for small (j) and in this way continue 

a 

the solution of (3) for high latitudes to low latitudes by solving 



2 ->->■ 

f^ V 4^ -f • Vf = gW (4) 

for small angles (J) , It is the purpose of this note to study the solutions 
of homogeneous versions of (1) and (4) in hopes of saying something about 
the solutions of the linear balance equation and methods for solving it. 
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We are interested in the homogeneous equation 

^{+ 

cos (j) 6X ) 



I ^ 

I cos (j) 6^ 



‘J’ I (cos (p ^ ) + ■ — J+ cos p ^ = 0 



£l 

si> 



( 5 ) 



(where we have dropped the 2Q/a ) all the way around the globe, so 
we may expect periodicity on X . We thus let 



4» = V(<p) 



( 6 ) 



We also will change the independent variable c|) by the obvious choice 

x = sin(|), (7) 

and write V((J)) E y(x). Then (5) becomes 



X 



I dx 








( 8 ) 



Suided by a suggestion in [1], we make one further change of variables 

, n/2 

y(x) = (1-x ) u(x) 



(9) 



to get 

2 

x(l-x^) ^ + [1 - (2n+3)x^] 4^ - (n^+2n)xu = 0 . (10) 

, Z QX 

dx 

Now (10) turns out to be an equation studied almost 100 years ago by 

Heun (see [1]), and is known as Heun’s equation, or occasionally as 

Lame’s equation. (For reference we list Babster’s parameters [1] for 

(10): a = -1, b = 0, a = n, 3 = n+2, y = 1, 5 = n+1, e = n+1 ). It is 

known that (10) has regular singular points at x = ± 1 and x = ^ and x = 0 

Using the method of Frobenius [2] we can generate modified power series 

about each of these points of the form 



(x - 






m=o 



c (x - 
m 



’‘o' 



m 



About the origin x — 0 (_the equator) the indicial equation for r has a 
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double root of zero. About the points x = ± 1 (the North and South 
poles) the indicial equation has roots r = 0, r = -n • 

The first zero value for r at x = 0 means that there is a solution 
of (10) about X = 0 which is an ordinary power series. This solution 
turns out to be an even function of x and is given by 



U = Uj^(x) 




C X 

m 



2m 



^nH-1 



(m + ^/2)(m + 1 + ^ ! 2) 

(m + 1)^ 



( 11 ) 

( 12 ) 



The ratio test applied to (12) shows that (11) converges for |x| < 1, 
and the test fails for |x| =1. A more subtle analysis shows that uj^ (x) 
converges for |x| = 1. The second zero value for r at x = 0 means 
that there is a logrcthmic solution near x = 0, which can be found by 
the method of reduction of order [2]. This solution is given by 



U 2 (x) = u^(x) 



/ 



’‘(-1)"'^^ dz 

2 2 
z(l-z^) u^(z) 



(13) 



An analysis of (13) shows that ^ 2 ^^^ does have a logarithmic singularity 
at X = 0> and 5 since u^(x) is finite at x = ± 1, ^ 2 ^^^ singular 

there, behaving like (1 ± x) ^ 

This is in agreement with the indicial equation for (10) at x = ± 1, 
which says that there is one solution which behaves as a constant (the r = 0) 
and one which is singular like (1 ± x) ^ . Thus the general solution 
to (8) can be written 

y(x) = (1 - x^) (c^u^(x) + C 2 U 2 (x)) (14) 
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where u^(x) is finite everywhere. Thus ^ 2 ^^^ singular at x = 0 

and , at x = + 1^ is more singular than (1-x^) . Thus the physically 

relevant version of (14) must have = 0» so that 

inX n /2 

® (cos (()) c u (sin (()) (15) 

n 1 

—00 

where the are the complex amplitudes of the longitudinal harmonics. 

Equation (15) is the solution to the homogeneous equation (5) and 
is not really useful, in itself, for solving the linear balance equation 
(1) . However, (14) tells us a great deat about solving (1) numerically. 

If one started at the poles or the equator (x = ± 1 or x = 0) with a 
known (finite) value of ^ then is necessarily zero. Then taking 

the exact solution, c^ remains zero. However, integrating numerically , 
one does not have the exact solution, and the roundoff errors effectively 
introduce a small amount of u^ into the solution. This would remain 
small if u^ were bounded. But integrating towards the equator (or the 
poles) the unphysical solution U 2 must eventually dominate. In other 
words, no matter how careful an integration scheme is used, because of finite 
precision in any computer, there is no way to integrate the balance equa- 
tion (1) from the poles to the equator and have a realistic solution at 
the equator . 

This has been instinctively recognized, a common solution method for 

the Northern Hemisphere models I 3] has been to stop integrating (1) at about 

20 degrees latitude and instead solve 

+ Vf • = g($) (16) 

O 

where k is sin(20 ). What can one say about the homogeneous version 
of (16)? Using the same substitutions as above we obtain, for the homo- 
geneous equation. 
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(l-x^)u" 



(17) 



- + 2k(n+l)x - ij u' - -^(x + kn + n)u = 0 



Equation (17) is not a standard equation, but is amenable to a standard 
power series expansion about the origin 

u (18) 

We then obtain the recursion relations 



6a^ ^ ” k ^2 ^1 ^ 

m(m-l)a = - — a ^ + fn^ + (2m-3)n + (m^-3m-2) a ^ 

m k m-1 L J ni-2 

k ra-3 

These involved equations show two things. First, if k is small (con- 
siderably less than 1), they tend to simplify, and this becomes a "singular 
perturbation" difference equation. Clhis is an almost unknown area which 
the author intends to investigate later.) Secondly, we note that no solu- 
tion of (17) is inherently odd or even, in contrast to both solutions of 

(18) which are even (see (11) and (13) ). Thus, although any function can 

\ 

always be written as the sum of an odd and an even function 



f(x) - J |f(x) + f(-x) I + "I |f(x) - f(-x)| , 

numerical integration of (16) will always introduce some odd elements of 
the solution. Since the equation one is really trying to solve is inherently 
even, these portions of the solution are unphysical. Since all solutions of 
(17) are bounded at x = 0 > the unwanted portion of the solution does not 
grow drasticly, so this is not as serious as trying to integrate (8). How- 
ever, it is expected that a singular perturbation analysis of the difference 
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equations (19) would show that the smaller the value of k , the more 
closely the solution to (17) will display the purely even character 
which (8) has. 

The non-homogenous versions (1) and (16) could be solved in power 
series in x = sin <p and Fourier series in X. The resultant double series 
show exactly the features above. With care, obviously an analytic solution 
in terms of a double series is possible. However, it certainly would be 
difficult to get reasonable answers from such an approach. The purpose of 
this exercise was to demonstrate that the exact solutions have inherent 
problems and thus any attempt to directly integrate the equations will give 
difficulty. 

The author would like to suggest an alternate approach to this problem 
which ought to be more stable, and which, in other applications, is often 
faster - the finite element method. This method has been applied to some 
elementary meteorological problems [4]. The basic idea is as follows: 
instead of taking an infinite series valid over the entire domain, one ex- 
pands ¥ as a collection of low order polynomials - a different polynomial 
for each grid rectangle in some discretization of the globe. The coefficients 
for such a ploynomial are determined by an integration of the trial solution, 
times suitable functions, over the entire globe. The justification for the 

computations is either an orthogonalization one (Galerkin) or a variational 
calculus one (Rayleigh-Ritz) . The latter is more in keeping with the use 
of NVA (numerical variational analysis) now used in meteorology, although it 
will cause a slight problem later. 
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For a given function $ it is easy to see that (1) is the Euler 



equation for 




+ 2g(4>)4< cos (J, d(|) dA . 



( 20 ) 



where D is the portion of the globe of interest. If we superimpose on D 
a set of grid points and connect these to cover D with triangles (as opposed 
to the rectangles normally used in finite difference solutions of (1)) we 
have a set of domains D. , In each of these domains we express ^ as a 



simple function (a low order polynomial) whose coefficients are to be 

determined. Some of the coefficients are determined by the need for con- 



The remaining coefficients are determined by minimizing (20) . The result 
is a matrix equation for these coefficients which is similar, but not 
identical, to the finite difference equations for (1). In particular, the 
fact that one is integrating (20) over an area means that the vanishing of 
f on a line (the equator) causes less problem than before. 

Experiments run with problems in stress-strain mechanics indicate that 
using linear functions for the H^i, or perhaps quadratics, will give excellent 
accuracy. The author has not yet run such an experiment for (20). 

A question with using (20) (ref erred to before) is that physically 
(1) may not be the best equation. to use from the poles to the equator. 
Equation (1) is actually 



This is really one equation on 2 dependent variables and 4> . In the high 

latitudes ^ is assumed known and one solves for ¥ , as in (1). In the 
equatorial latitudes however, it is believed [5] that ¥ is better known and 



1 



tinuity of ¥ and some of its derivatives across the boundaries of D^. 



2 2 

fV f + Vf • V4' = V $ 



( 21 ) 
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one uses C21) to solve for 'J> . TKe variational formulation of C21) 

for $ is 






-If 



[(vl)^ + 2h ('!')$-] a^cos 



C22) 



where h(4^) is the left side of (21). Then one could take tlie domain for 
(20) to be down to some reasonable latitude, say 20 degrees N, and then 

o 

use (22) through the equator to 20 S., then use (20) down to the South 
pole. The difficulty is that ^ must be precisely specified on the boundary 
(20 ) for (20) and <l> there likewise for (22). That is not very realistic. 

Haltiner has proposed [5] an iterative method to avoid the effects of 
inaccuracies in specifying the boundary conditions. It basicly involves 
solving three problems in overlapping regions. It involves solving for ^ 

o o 

from, say, 40 N. to 40 S., using some reasonable approximation (say ^/f) 

for on the boundaries. This computation is done by using the observed 

2 

winds to find the vorticity C and solving V = C • Tke resultant values 

o o 

of at, say , 20 are used to solve (1) for from the poles to 20 

o o 

Using these values of ^ from, say, 30 to 20 and the values of ^ from 

O o 0.0 

±20 to 0 from step one, one can then solve for $ from 30 N. to 30 S. 
using the linear (or non-linear) balance equation. Because of the inherent 
character of the Laplacean as a smoothing operator, this overlapping pro- 
cedure should give more reasonable answers. 

Another approach which minimizes this dichotomy is to use a finite 
element-Galerkin approach, a non-variational method. This requires that we 
use (21), the differential equation, rather than (20) or (22). One then 
subdivides the entire globe into and takes $ and ^ to be low order 
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polynomials in each . In those domains where ^ is believed known, 
the coefficients of are specified. Then the coefficient of the re- 
maining functions are determined by insisting that, for all i , 




1 



v^. 



Vf - 



v^$, 



p. dA = 0 
1 



( 23 ) 



1 

where the p^ are the low order basis polynomials. This is nearly 
equivalent to the previous two integrals, but conceptially places both 
functions 'F and $ on the same basis. The result of (23) is again a 
linear matrix system for the coefficients of the and , the poly- 

nomial functions which represent ^ and 'F in each grid area. 

A potential advantage to this approach is that there is no need what- 
soever to make the grid spacing rectangular or uniform, because there is no 

finite differences whatsoever. The and are explicit functions, 

11 

and are differentiated exactly by hand. By coverting to the x , X 

coordinate system (x = sin <()) , the coefficients in (23) are polynomials, 

so that the computations and integrations are elementary and in fact are 

exact. Thus the computational errors are only in the use of the and 

'F. , and in the data itself. 

1 
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